Quantum Monte Carlo

Video of SciPy 2011 talk

2011-10-05T23:18:00.000-05:00

The video from my SciPy 2011 talk is available here.

At the end of the talk, Andy Terrell asked a question about why this structured approach is necessary - why not just program the steps directly in a script?

I agree that the system I presented does seem very 'heavyweight'. But it feels that doing the transformations directly loses some meta-information that would be useful to keep around, such as which transformations were applied at each step. Additionally, by storing the the pieces of the derivation in data structures, they are potentially easier to manipulate (and emphasizes that they can and should be programmatically manipulated). However, Python's meta-tools are quite powerful, and could probably used to manipulate derivations written directly in a script.

Switching to a different issue, the presented example effectively inlines all the functions into the integral. This is clearly not scalable in the face of larger and more complex expressions. The system needs some way to represent pieces that should transform into a function call in the generated code. I have something working now, but it needs some clean-up work.

Modeling derivations: Questions and Answers

2011-06-07T07:32:00.003-05:00

The past two posts (part I and part II) described a system (or style) for scientific programming that I've been working out. This post will have more in a Q&A format (The questions are mostly ones I ask myself - perhaps it should better be called a 'rhetoric & rant' format) An abstract for a talk and paper was accepted to SciPy 2011, so hopefully I can explain it clearly enough by then to get some feedback.

What is the target application? Goals?

The target goals are better enabling algorithm research and development (that is, writing a scientific code), and a flexibility to target different runtime systems.

Looking at the examples and case studies on various computer algebra system (CAS) websites, I see many cases where the user modeled some system, solved it, and then used the solution to improve or fix the system. The hard part seems to be modeling the system - this is not my target. I'm interested in problems where the second step is the hard part - and one needs to research new and better solvers.
I'm thinking in terms of a traditional scientific code, and how to automate parts of the development workflow. One big issue is how to use a CAS to assist in the development. Also I'm interested in problems that require high performance computing (HPC) - basically problems that will consume as much computing power and algorithm research as possible. This likely means the code will need to be tuned and altered to support different target systems (multicore workstations, distributed memory parallel clusters, accelerators like GPU's, etc).
Why not use a compute algebra system directly? It looks like you're creating a lot of extra work for yourself?
I would like a separation between the specification of an algorithm (for scientific computation, this is largely expressed in equations) and its implementation in the runtime system (along with an automated transformation from the spec to the final system).

Typically a CAS is associated with full featured programming language, but implementing the algorithm in that language doesn't preserve this separation.

As a concrete example, a CAS typically includes facilities for performing numerical integration. What if I want to develop a new numerical integration method - are there any examples on how a CAS can be used to do that?

Another issue with implementing HPC-style problems in a CAS is what happens when the existing solvers in the CAS are too slow, or the CAS system language is too slow (eg, how to scale to large parallel systems). The CAS may do a great deal of work for you, but what happens when you want to move past it? This is the part where code generation comes in.
What is the role of code generation in this system?
More generally, the 'code generation' step is the 'transformation to the target runtime system' step. The transformation may be trivial, if the final equations can be solved by the CAS facilities, but I'm anticipating the need to ultimately target C/C++ or Fortran.
Also, it seems useful to introduce this method in small pieces into an existing code (converting an entire code would usually be too dramatic of a change to be successful). It would smooth the transition for the generated code to fit the style of the existing code. A template system for generation would be a useful approach. A text-oriented (or basic token oriented system, like the C pre-processor) would work, but they have problems. A system oriented around an actual parse tree of the target language would allow the most powerful transformations.

Is this related to reproducible research in any way?
Packages for reproducible research seem to be focused on the level of capturing parameters that are input to a scientific code, so that the results can be repeated and revisited later.
Executing a program may be reproducible, but the construction of that program is not.

Aside:
Reproducibility can be achieved through documentation or automation. The previous standard for scientific work depended primarily on documentation. With computers, automation becomes possible as a working method. (For space reasons, this analysis is very simplistic --- there is more to it.)
</aside>

This work is intended to make the construction of a code automated. One advantage is in making corrections to mistakes in the derivation - fix it, push a button, and the final code is automatically updated.

This should be a complementary to other systems for reproducible research.
Why are you doing this?
I'm tired of tracking down missing minus signs and misplaced factors of 2 in my code. A significant feature of this approach is that the entire chain of derivation and generation of the program is performed by machine, to eliminate transcription errors. (If this sounds too rosy, rest assured there are plenty of other types of errors to make).

Example of modeling derivations

2011-04-23T23:19:00.002-05:00

The last post discussed a style of scientific programming that modeled derivations symbolically to create a machine-checked chain from starting equation to final code.

This post will provide an example using the configuration integral for the canonical ensemble in statistical mechanics. I want to evaluate the integral using various methods, starting with grid-based integration methods. Using a grid-based method suffers from the curse of dimensionality, where adding particles raises the time to compute exponentially (which is why Monte Carlo methods are normally used). However, it can still be useful to check Monte Carlo codes with other integration methods for small particle numbers. Also, the free energy (and all associated quantities) is computed directly, where it is more difficult to compute with MC methods. Finally, I'm curious as to how large of a system could be computed on today's hardware.

Sympy is the computer algebra system I'm using. The code for this example is available on github in the derivation_modeling branch, in the prototype subdirectory.

The output (MathML and generated python) can be seen here.

The flow is as follows:

Derive the composite trapezoidal rule, starting from the expression for a single interval

Generate code for the composite trapezoidal rule

Start from the configuration integral (which I've called the partition function in the documents) and transform it by specializing various parameters - 2 particles in 2 spatial dimensions, using the Lennard-Jones interaction potential, etc - until it becomes completely specified numerically.

Generate code that evaluates this integral using the generated integration code from step 2

There's still a tremendous amount of work to do on this, but hopefully an outline of this style of programming is visible from the example.

The MathML pages have a lot of detailed small steps - it would be nice to have a collapsible hierarchy so all the little steps can be hidden.

Some abstractions should remain during the steps - in particular the potential. The expression for the integrand can get quite messy from explicitly inserting the expression for the potential (and this is a simple case - more complex potentials would make it unreadable).

More backends for the code generation - C/C++ is the obvious choice in pursuit of higher performance. (The generated python can be run through Shed Skin, for another route to C++). Converting to Theano expressions might be a quick way to try it out on the GPU.

Modeling derivations

2011-04-07T01:55:00.002-05:00

A previous post describes how I kept reference versions of important derivations in a document. I would like to move the idea out of LaTeX and into some format more amenable to machine checking. To do this requires some structure for describing and modeling a derivation.

A derivation starts with some base equation, then proceeds through a series of steps transforming it into a form that can be solved numerically (A symbolic/analytic solution would be nice, but not possible in most cases.) The steps can be divided into three categories

Exact transformations - rearranging terms, replacing an expression with a different form, the usual steps involved in a proof, etc

Approximations - using a truncated series expansion, replacing a derivative with discretized version, etc.

Specializations - Specifying the physical or model parameters in order to make the equation numerically solvable - like the number of spatial dimensions, number of particles, interaction potentials, etc.

At the end of the derivation, the result can be evaluated numerically using the CAS (Computer Algebra System - I assume this derivation modeling has been implemented using one) facilities, if simple enough. More complex or time-consuming results (a more likely case) can be further transformed by code generation to C/Fortran (or some other execution system).

One important feature of this workflow is every step in creating the program is checked by machine. This should greatly reduce transcription errors (dropped minus signs and factors of 2 - this is what started the whole project in the first place). Also it moves part of the process of creating scientific programs into the realm of 'ordinary software', and as such makes it subject to normal software engineering techniques - unit testing, etc. Testing scientific (numerical) software is a difficult problem - this should help make it (a little) easier and more reliable.

Some other work I've run across related to describing mathematics:

Structured derivations were developed for teaching proof in high school mathematics.

OMDoc is a markup language for describing mathematical structures at a higher level than a single formula (a proof, for instance)

Levy constrained search in VMC?

2009-08-27T00:19:00.000-05:00

The Levy constrained search formalism is a useful way to reason about density functional theory. It finds the minimum energy using two nested searches. The inner search seeks the minimum energy wavefunction for a particular density. The outer search seeks the density with the minimum energy.

This could be directly implemented using VMC. The main ingredient is energy minimization subject to a density constraint. Has anyone done this?

It wouldn't be very efficient compared to optimizing the VMC energy directly, but it would make an interesting connection between the methods.

Why scientific computations fail

2008-12-14T14:49:00.002-06:00

(I was busy this fall with another project. Now it's time to return to scientific programming)

The previous post covered several levels of tests for assessing the correctness of scientific computation. Now let us look at reasons why the results might not match what is expected.

Programming errors
Mistakes in derivation and transcription of equations
Violations of assumptions or conventions
Errors from floating point
Errors and undesirable behavior from the algorithm (often problem-dependent)
- Algorithms that have been well-studied by numerical analysts
- Those that have not

Programming errors

The usual sorts of problems in programs (array indices, memory problems, etc). This category of problems can be tested with unit tests and can make use of tools and practices from mainstream software engineering.

Mistakes in derivation and transcription of equations

This is similar to the implementation not matching the specification in traditional software engineering. With scientific computation, though, the "specification" could be more precise in principle (equations and algorithms), and computer algebra systems could aid with performing and verifying steps in the derivation.

Violations of assumptions or conventions

Usually happens during derivation of the equations. A stronger connection between the derivation and the final code might help this.

Errors and instabilities from floating point

Accumulated round-off, loss of precision from subtracting nearly equal numbers, and other problems. For more subtleties than you every wanted to know about, visit W. Kahan's site.
One of the most annoying parts about floating point is lack of good model for understanding (and analyzing) floating point errors. I will return to this issue in a future post.

Errors and instabilities from the algorithm.

This category might be better divided into two subcategories - algorithms that have been well-studied by numerical analysts (like matrix operations, certain ODE's and PDE's, optimization), and those that have not. Many algorithms have been well-studied and have nice packages available to use as black-boxes. However, most scientific programs are a composition of well-known algorithms along with some problem-dependent features as well.

Often times the behavior of the algorithm depends on the details of the physical system (ie, is it a metal or insulator)?

This category also includes science-based and numerical approximations.

Now, how do we go from test case failure to deciding which of these categories the problem falls into (and then, of course, to locating the precise problem)?

Unit testing will help with the first category. But not so much with the rest. Unit testing is helpful when you want to make sure things haven't changed. The hard part is getting it right the first time.
The different levels of tests in the previous post should give some insight into which category the problem might occur.

First I want to look at various models of floating point error and see if we can gain some understanding of the intrinsic uncertainty in a calculation. This seems to be a necessary prerequisite for answering the question of 'how close is close enough?', in order to decide whether a test has actually failed or not.

Assessing the correctness of scientific computation

2008-07-13T21:54:00.000-05:00

In the previous post I promised some thoughts on assessing the correctness of scientific computation.

Computational science programs are typically assessed with a series of test cases, ranging from order-of-magnitude estimates ("those numbers don't look right") to comparing with experimental data. I've attempted to divide these test cases into categories, where each category shares similar expectations as to what 'close enough' means and has similar implications into potential causes for discrepancies.

Order of magnitude estimates from general principles, back of the envelope calculations, or experience with similar problems.
Analytically solvable cases.
Small cases that are accurately solvable
For small systems, or simple cases, there are often multiple solution methods, and some of them are very accurate (or exact). The methods can be compared for these small systems. An important feature of this category is any approximations can be controlled and the effects made arbitrarily small.
Results from the same (or similar) algorithm
Comparing the same system with the same algorithm should yield close results, but now there is additional uncertainty in the exact implementation, and the sensitivity to input precision (particularly when comparing results from a journal article)
Results from a different algorithm
This looks similar to the second category (comparing with exact or more precise algorithms), but this is the case where there a multiple methods (with different approximations) for handling larger systems. There are likely to be more approximations (and probably uncontrolled approximations) involved.

(In electronic structure there is Hartree-Fock(HF), several post-HF methods, Density Functional Theory, and Quantum Monte Carlo. And possibly a few others.) Now one has to deal with the all the above possibilities for two programs, rather than just one.
Experimental results
In some ways this is quite similar to the the previous category, except determining "sufficiently" close requires some understanding of the experimental uncertainties and corrections in addition to possible program errors.

The testing process for each case involves running the test case, deciding whether the program result is 'close enough'. If it is, proceed to the next test case. If it is not, work to discover the cause of the discrepancy. Next post I'll look at a list of possible causes.

Engineering scientific software

2008-07-12T01:20:00.002-05:00

Some pointers from Greg Wilson on engineering scientific software.

The first pointer is to the "First International Workshop on Software Engineering for Computational Science and Engineering".

A couple of Wilson's short summaries of papers from this workshop
(from recent research reading)

Diane Kelly and Rebecca Sanders: “Assessing the Quality of Scientific Software” (SE-CSE workshop, 2008). A short, useful summary of scientists’ quality assurance practices. As many people have noted, most commercial tools aren’t helpful when (a) cosmetic rearrangement of the code changes the output, and (b) the whole reason you’re writing a program is that you can’t figure out what the answer ought to be.

Judith Segal: “Models of Scientific Software Development”. Segal has spent the last few years doing field studies of scientists building software, and based on those has developed a model of their processes that is distinct in several ways from both classical and agile models.

What I found interesting:

Kelly and Sanders suggest that the split between verification and validation is not useful, and propose 'assessment' be used as an umbrella term. (I will return to the point in a later post)

The second pointer is to the presentations at TACC Scientific Software Days.

One of the presentations there was by Greg Wilson himself on "HPC Considered Harmful". There is one slide on testing floating point code (slide 15). In particular the bullet point testing against a tolerance is "superstition, not science". (And he's raised this point elsewhere)

This statement has bothered me since I first read it, and so the next few posts will explore some of context surrounding this question. The first post will list different methods for assessing the correctness of a program, the next post will categorize possible errors, and finally a post looking at 'close enough'.

Computing an unbiased inverse

2008-05-29T22:23:00.000-05:00

In a previous entry, Reweighting is biased, I had discussed the problem of a finite-sample size bias in reweighting. In this entry I will investigate a possible solution.

But first, I had asked if this was known in the literature. A comment by Lee Warren pointed to this article, Population size bias in descendant-weighted diffusion quantum Monte Carlo simulations, which discusses this very issue (in a slightly different context, but the root issue is the same). In that context the problem is even more severe since the straightforward solution (increase N) carries a heavier computational cost.

One of the crucial issues is finding a good estimator for the inverse (1/D). The naive estimator (simply do the division) is biased. So the question is: are there any other ways of performing division that could be used? Yes, in digital arithmetic there are a couple of algorithms for division using only multiplication and addition - Newton-Raphson and Goldschmidt (see the Wikipedia entry on digital division)

In general, to have an unbiased estimator, each sample of the random variable must enter the estimator linearly (alternately, we can say each sample should only be used once and then discarded)

The Newton-Raphson iteration for the inverse (1/D) is x_i+1 = x_i(2-D x_i). To make this unbiased, in each iteration the value of D needs to be an independent sample, and the values of x_i from previous iterations also need to be computed from independent samples. One can think of the algorithm forming a tree, with each iteration branching to two children (each child is one of the x_i computed by the previous iteration), until reaching the leaves representing the initial iteration (x₀ = T₁ + T₂ D)

I tested this algorithm, and it appears to give unbiased results. There is a catch, though. The result has a large variance. Using the average of all the samples with the naive estimator yields a bias that is much smaller than the variance of the unbiased estimator (hence the mean squared error is a more appropriate metric for the combination of variance and bias)

Open questions/issues:

Is there some modification of the algorithm that will produce a smaller MSE than the naive estimator? Or even an MSE that isn't too much larger?

How many iterations are necessary - how to analyze the convergence of this stochastic Newton-Raphson iteration?

Try out the Goldschmidt algorithm.

Reproducible Research

2008-04-06T00:56:00.000-05:00

There were a couple of recent posts on Greg Wilsons's blog about reproducible research. The first is a pointer to a paper about WaveLab. The paper is from 1995, although it appears the philosophy hasn't changed much looking at a more recent (2005) introduction on the website. The second post contains a reference to the Madagascar project.

These projects follow the work of Jon Claerbout on Reproducible Research.

Thoughts

This work is all about designing a workflow for computational research, and creating tools to support that workflow. Much of it seems quite banal, but part of the whole purpose of doing this is to free you from having to think about the boring and repetitious steps.
Levels of reproducibility
1. Yourself, right away (when tweaking figures, for instance)
2. Another research group member, or yourself 6 months later
3. Someone in another group, (a) able to run the code and get the same results, and (b) able to understand the code
4. Someone in another group, able to use your description to write their own code that gets the same results
The last level is the gold standard in scientific reproducibility. Experiments are reproduced from published descriptions (ideally, although the descriptions aren't always complete). Or theorists rederive intermediate steps leading to a paper's conclusion.

In software development, levels 1-3 are what counts. (It's not very useful software development if others can't run your code and get expected results.). The last level is similar to multiple implementations of a standard (eg, C++ compilers from multiple vendors.)

The referenced works seem to focus mostly on processing of raw data (images, seismic imaging, etc) and the processing steps to produce a final result.

How would this apply in other fields? Say, ones with smaller input data sets, but much higher computational demands.

The Madagascar project, for example, uses SCons to coordinate the steps in producing research output. What happens when one of the steps involves several weeks of supercomputer time? Two immediate problems that spring to mind - one probably doesn't want this step launched without warning by the script. Furthermore, the steps for launching such a job are likely to be very site-specific, making it difficult for others to reproduce elsewhere.

Somewhat related entries of mine: Keeping notes and Solving Science Problems

Notes from PyCon 2008

2008-03-19T22:50:00.001-05:00

I attended PyCon 2008 this past weekend, and had a good time meeting people in the Python commmunity.

Some notes

I tried to give a lightning (5 minute) talk about sympy and code generation. Unfortunately, Ubuntu 7.10 on my laptop (Dell Vostro) failed to display to the external output, so I was unable to present to the larger audience. However, I was able to give the talk to the people at the science BOF session. The slides are available here

At the science BOF session I met Michael Tobis, the author an interesting system called PyNSol. It uses python to describe the differential equations, and generates Fortran for performance.

The python 2 to 3 converter looks interesting. Not so much for the conversion of python versions, but for the infrastructure developed for the tool.

Enthought was showing off some multi-touch displays. Being able to interact with 3D visualizations (and Google Earth) via gestures for zoom, pan, and rotation made for a nice demo. The real question is if this type of display would be useful for actual use.

Code generation now uses templates

2008-03-17T20:34:00.000-05:00

In Quameon, the derivatives for the Gaussian basis functions are computed with sympy and the python code is generated. Creating the code that surrounds the generated code is a bit of a tedious procedure - the syntax tree needs to be constructed by hand. But no more! The surrounding code can now written in python and a parser (written using pyparsing) will create the syntax tree for you.

Some details:

Identifiers surrounded by percent signs (eg %name%) get turned into a template node type for subsequent replacement.

Example template:


class %orb_name%:
  def compute_value(self,alpha,v,r2):
    x=v[0]
    y=v[1]
    z=v[2]
    %value%
    return val

The code replaces %orb_name% with the orbital name (s,px,py,px,etc) and %value% gets replaced with the expression for the value of the orbital (the methods for the derivatives not shown - see codegen/primitive_gaussian/deriv.py in SVN for the full template.)

Why python?

2007-12-21T08:31:00.000-06:00

The author of the first comment on the Quameon entry wonders if Python will be too slow.

One reason to use Python is ease of programming. Hopefully the advantages of programming in Python will outweigh the performance penalty.

However, at some point, runtime speed will be an issue. Here's a list of steps to improve runtime performance (in increasing order of effort/expected return)

Use psyco. Put "import psyco; pysco.full()" at the beginning of your program and get an instant speed increase.
Use the Shed-Skin Python to C++ compiler. (I haven't tried this yet - don't know if it works with Quameon.)
Profile and rewrite slow parts in C/C++/Fortran. I've done this with the inter-particle distance computation in a classical MC code in Python, and gained some performance.
Code generation. Right now, the code for the Gaussian basis functions is auto-generated (no more computing derivatives by hand!!). I would like to move more of the basic formulas to this technique. Then it should be easier to retarget the code generation to another language (in theory, I haven't tried it yet).

Other tidbits:

Monte Python uses Python and C++ for QMC (paper, code). They wrote the time-consuming parts in C++, and use Python to glue it all together. The algorithm variants they tested (parallel distribution strategies for DMC) were achieved through changing the Python part of the code.

Code generation could potentially result in *faster* code than writing a general framework in C++ or Fortran. This is because the code can be specialized for a given physical system, exposing more optimization opportunities. For an example of code generation, look at the development SVN area of PyQuante - there is a program that will produce a C++ program to compute the energy for a molecule given in the input file.

Quameon

2007-12-05T23:32:00.001-06:00

I've been working on a Quantum Monte Carlo code implemented in Python. It's called "Quameon" (the name is a mixture of letters from "Quantum Monte Carlo in Python"). There's a SourceForge project for it - here are links to the SourceForge project page, the project web page, and the project wiki.

The code doesn't do much currently, but I mention the project now so I can write blog entries about various aspects of getting the code working.

I'm currently trying to validate that the code works by using HF single particle orbitals (obtained from PyQuante) and no jastrow factor - this should reproduce the HF energy. It seems to work for a few atoms (He, Be, C), after sorting out some normalization issues with the basis sets (apparently all basis sets are not normalized - this doesn't affect the energy, but does affect the eigenvector coefficients). I will need longer runs to be sure. The next step is to profile the code to look for any quick and easy performance tuning opportunities.

The first target (after getting the code running reliably) is to implement several forms for jastrow factors and several VMC parameter optimization algorithms.

Keeping notes

2007-10-04T22:38:00.000-05:00

Keeping a record of work is important for research. I use two main tools for this purpose. First, a text file with entries in chronological order (one file per month, with names like "Sept2007"). I write the date and then start taking notes for that day in the file.

Pro - easy to use, chronological view easily available
Con - math hard to write well, ideas spread over several days hard to view topically.

The other tool I use is a wiki. It's useful to keep data in a structured, hierarchical format that's easy to edit.

Pro - ideas arranged topically, can do better looking math
Con - must work topically from the start, no good chronological view (Yes, I know you can view changes in order, but it's not as easy to reconstruct what I was thinking on a particular day.)

I want a single tool that does both of these - that can view the notes both chronologically and by topic.
For work-flow, I want the ability to enter ideas and notes free-form, and perform categorization and annotation at a separate time.

Another additional feature that would be nice is recording other actions on the system (at least by reference).
One example is check-ins to source control.

Another possible example - What if gnuplot were to record all the commands issued (and kept a copy of all the files used?) Graphs could be re-created and played back at a later date easily.

I looked around a little bit to see if there are other applications that might meet some of these requirements (or at least be an improvement over text files and a wiki). One class of app is a note taking application. I looked at Tomboy. It does have a Latex plugin (though I could not get it to work). Latex side-rant: It's great for making nice looking mathematics, but you wind up with "equations under glass". It's pretty to look at, but you can't manipulate it further. For mathematical content, I want something with more precise semantics (Content MathML, for example)

Another class of tool is mind mapping software. I looked at Vym (View your mind) briefly, but I don't think I will like this type of software. Viewing a graph of links is nice, but not as the main working view - I prefer a wiki view (pages with links)

Black Box Reweighting

2007-08-19T21:56:00.000-05:00

In the last post, I wondered whether the black box reweighting method might help with the bias.

The Black Box Reweighting (BBRW) method computes weights using the observed distribution of samples, rather than the targeted distribution. (The author's website contains a link to supplementary material with further discussion.)

I ran a simple test, and the bias was slightly worse with the BBRW method, compared to standard reweighting.

The BBRW method does a fine job of solving the problem the authors intended (correct handling of a distribution that may be incompletely sampled, and lower variance of the result - the lower variance of the result was apparent in my tests), but it doesn't solve the reweighting bias issue.

The search for a method to correct reweighting bias continues...

Reweighting is biased

2007-08-14T22:52:00.000-05:00

The reweighting technique, where samples from one probability distribution are adjusted with a weight factor to compute averages from a different (but usually very similar) distribution, is biased for a finite number of samples.

The details can be found here.

Question for the audience: Is this known? I don't recall seeing this mentioned any time the method has been presented.

Some possible routes to fix the problem are listed in the paper, but alas, I haven't been able to turn any of them into a useful solution yet.

VMC optimization papers

2007-08-10T22:29:00.000-05:00

A recent comment suggested these papers on VMC optimization:

Alleviation of the Fermion-sign problem by optimization of many-body wave functions

Wave function optimization in the variational Monte Carlo method

I'd also add these to the pile of papers to read and understand:

Optimization of quantum Monte Carlo wave functions by energy minimization

A simple and efficient approach to the optimization of correlated wave functions

And since this post is alread about VMC optimization, some further questions:

What about simultaneous geometry and parameter optimization?

The focus seems to be on single geometries - what's likely more interesting in the future is a family of geometries - has anyone parameterized the VMC parameters (by the bond lengths or nuclear coordinates)?

Correct estimators in DMC without forward walking?

2007-08-10T21:09:00.000-05:00

Computing properties other than the energy using Diffusion Monte Carlo requires forward walking to get correct answers. That may change, if this paper - Hellman-Feynman operator sampling in Diffusion Monte Carlo calculations - is correct.

Although as I look at the paper more, the technique doesn't look any easier than the forward walking implementation by Casulleras and Boronat (the paper does mention this similarity). It appears to use similar data in a different combination. I would be interested to see how the noise, bias, and stability characteristics compare.

The Hellman-Feynman theorem seems to be getting lots of attention in the QMC world recently. It's also been used to construct Zero-Variance Zero-Bias estimators (most recently example here).

Gravitational N-body and the Art of Computational Science

2007-08-02T22:35:00.000-05:00

A rather ambitious project to describe a gravitional N-body code and how to write it in Ruby: The Art of Computational Science

As Greg Wilson points out, something like this would be nice to have in Python.

Efficient Global Optimization

2007-03-17T23:51:00.000-05:00

The optimization method mentioned last post is an improvement on Efficient Global Optimization, presented in the paper Efficient Global Optimization of ExpensiveBlack-Box Functions (also available from on one of the author's websites)

Another paper listed in the references, A Taxonomy of Global Optimization Methods Based on Response Surfaces(by D. R. Jones, one of the authors on the EGO paper) is a nice overview of various optimization methods. It presents each method and gives ways in which it can be fooled into choosing the wrong point as the optimum. A fix is given, which then leads to the next method in the list.

One of the core ideas is regression based on Gaussian processes (kriging). To learn more about Gaussian processes, see www.gaussianprocess.org. The two Jones papers listed above also show how the GP based methods fit in with other basis set methods.

Useful paper for VMC optimization?

2007-03-02T17:44:00.000-06:00

I ran across this paper, An informational approach to the global optimization of expensive-to-evaluate functions, and wondered if it might be useful for VMC optimization?

The type of problem being solved seems to fit:

Function is expensive to evaluate.
Function evaluation may be noisy.

Many optimization methods estimate an optimal point and evaluate the function at that point. The method described in the paper estimates the uncertainty in the knowledge of the function, and evaluates the function where it will best improve our knowledge of the function.

Some potential issues:

This method doesn't make use of gradient information, which may make it less competitive than methods that do use gradient information.

How well does it scale with the number of parameters? The authors present a method to keep the cost under control as the dimension increases, but is the method still effective then?

March meeting

2007-02-23T08:21:00.000-06:00

This year's APS March meeting has a couple of sessions on QMC: H4: Recent Advances in quantum Monte Carlo Simulations and V21: General Theory: Computational Quantum Monte Carlo Methods.

Another session of interest (for me, anyway) is A23: Focus Session: High Pressure I - Earth and Planetary Materials

I will not be attending the March meeting this year, as my wife and I are expecting our first child in mid-March.

Science Blogging Conference

2007-01-18T08:24:00.000-06:00

The North Caroline Science Blogging Conference starts this weekend. It looks like it covers some introductory how-to material, and then dives into issues such as promoting public understanding of science, teaching, and how does blogging interact with research.

If you're put off by the term 'blog' being overused and overhyped, think of the conference as discussing the potential and challenges of a (relatively) new one-to-many and many-to-many communication mechanism and how that could affect the interactions between scientists and the public, and the interactions between scientists.

Basis sets galore

2007-01-17T21:52:00.000-06:00

More basis sets that you can shake a stick at over at the EMSL Basis Set Exchange. (Not sure why you would want to shake a stick at basis sets, though.)

Schrodinger equation in the mass media

2006-10-26T08:08:00.000-05:00

The Schrodinger equation appears in the background in the second half of Weird Al's video White & Nerdy.

The video also shows rolling dice (briefly), if one wishes to stretch its relevance to include Monte Carlo as well.

The Other QMC

2006-10-20T21:12:00.000-05:00

Quasi-Monte Carlo is the other method commonly referred to by the QMC acronym (I will abbreviate it QuasiMC to minimize confusion). It's not really Monte Carlo since the sequence of points is not random, but the point sets do have a Monte Carlo-ish feel about them.

And in fact, it does share the property of arbitrary termination with Monte Carlo integration (consider a grid - you have have use all the points at a particular spacing. Stopping half way across the region of integration is not going to yield good results).

The reason people are interested in QuasiMC is that its convergence is better than the 1/sqrt(N) of Monte Carlo. It does this by using point sets that are more evenly distributed than random points - there are fewer clumps of points or large gaps that appear in random sequences (Search for 'Low-discrepancy sequences'. Here's the Wikipedia entries for a technical discussion and some sequences).

QuasiMC point sets are usually generated from sequences where the low-discrepancy condition can be verified theoretically, but it is possible to use the intuitive approach of making the points maximally spread out (see this presentation for an example).

The following algorithm will generate a set of N points that gives better convergence than random (at least it worked on a simple integrand in one dimension.)

compute some random points (>> N points)

pick a starting point and put it in the set S.

find the point that is the furthest from all existing points in S, and place that point in S.

repeat the step 3 until N points are selected.

(I'm not saying this an efficient way to generate the points, just that it can be done.)

For addition information, see also chapter 7.7 in Numerical Recipes

QMC wiki

2006-09-01T20:21:00.000-05:00

Check out the QMC wiki at http://www.qmcwiki.org

It looks like a promising resource for the QMC research community.

Optimal histogram bin width

2006-06-27T01:15:00.000-05:00

Kevin Knuth wrote a paper about finding the optimal number of bins to represent data in a histogram (Optimal Data-Based Binning for Histograms). He starts from a piecewise constant density model and finds the (Bayesian) posterior probability from this model (equation 36, which is actually the log of the posterior). The posterior function is then maximized to find the number of bins that best models the data.

The article also investigates the number of data points for a reliable estimation of the density. The recommendation is 100-150 points, if the distribution is Gaussian.

It would be interesting to apply this method to radial distribution functions. However the assumption of a constant volume for each bin is not met in this case. There are several ways this could be adjusted, but I'm not sure they are valid (scale each bin count by the volume, or use non-uniform bin spacing to maintain constant volume)

Alternately, the discussion references other algorithms for dealing with variable bin-width models (which may be better for resolving multiple peaks anyway).

QMC derivation notes

2006-06-08T23:06:00.000-05:00

I posted a document I wrote in grad school, Notes on the wavefunction and local energy. It contains derivations of various QMC formulas, particularly the first and second derivatives for several forms of wavefunctions.

I'm posting this for two reasons. The first is in case anyone finds the formulas useful when working on a QMC code.

The second is related to the process of scientific programming. When writing a QMC code, I found it useful to record the formulas and derivations in a neatly typset form. Then the next step involved turning the equations into computer code. (Then, of course, testing and debugging).

This workflow is what I would like to capture with the Progamming in Mathematical Notation work. The document with derivations could be written in content MathML (or something more amenable to human manipulation). Ideally the computer could then assist with verifying the derivations for correctness, and with converting the equations into computer code.

And as long as I'm dreaming, I'd really like a wiki-like interface for creating and editing such a document (making a set of hyperlinked pages rather than a single linear document)

Interesting article studying writing journal articles

2006-04-05T22:43:00.000-05:00

Last fall I looked at methods for improving the efficiency of MC simulations, and was pointed to the article, Efficient Monte Carlo methods for the computer simulation of biological molecules by Djamal Bouzida, Shankar Kumar, and Robert Swendsen.

A search of the author's names on Yahoo turned up a fascinating article by Ann Beakeslee, who studied the process of writing this paper. The article examines the interaction of the graduate student (Bouzida) and the advisor (Swendsen), covers the process of learning to write scientific papers, and looks at difficulties Bouzida faced as a novice.

Search for chemicals

2006-03-22T23:24:00.000-06:00

Check out chmoogle, the search engine for chemicals. Many of the results are for suppliers (not so useful for theoreticians), but some of the results contain more information about the chemical.

The "details" link by the chemical structure picture also has alternate names, along with links to search Google or Yahoo.

[APS] First Principles Molecular Dynamics on Blue Gene

2006-03-20T23:13:00.000-06:00

Francois Gygi gave a talk about scaling a first principles molecular dynamics code to 65,536 processors on the Blue Gene/L at LLNL.
The core numerical algorithms that need to scale are the FFT and dense linear algebra. The FFT is limited in scalability to 512 processors per k-point. Fortunately, there are many k-points and they can be computed independently.

Modifying the assignment of tasks to processors increased the performance by 64%! FLOPS are free, it's the communication patterns that are determine performance.

One research problem is improving the scalability of computations for small systems (ie, 32 water molecules), so they can be simulated for longer times.

[APS] Solar Power

2006-03-20T22:46:00.000-06:00

David Carlson (from BP Solar) gave a talk titled, The Status and Outlook for the Photovoltaics Industry. His presentation slides are also available online.

Some random facts I found interesting

The photovoltaic industry has been growing at 35% per year.

This year, the PV industry will use as much or more silicon than the semiconductor industry. This may cause problems with the growth curve, as suppliers ramp up to supply the PV industry.

The current cost per kWh is $0.18. Less than half of the cost is due to the PV module. (Slide 30)

Given the cost of PV electricity, the cost of conventional electricity, and the availability of sun, Spain is almost to the point of cost-competetive PV electricity. (Slide 31)

[APS] QC for QC

2006-03-13T17:35:00.000-06:00

Intersting point from Alan's talk about using Quantum Computing to solve Quantum Chemistry problems: computing cholestrol is about as hard as breaking encryption algorithms.

The quantum computing approach takes quite a bit of setup - an actual Hartree-Fock calculation or more. The advantage of QC is that the exact answer (Full CI) can be obtained easily[1] by the QC from that starting point. (But going from HF to FCI is the hard part classically - I think Full CI scales exponentially in the system size)

[1] as much as anything with quantum computers can be called _easy_

[APS] Hydrogen storage

2006-03-13T16:47:00.000-06:00

I went to a few talks in a session about Simulations of hydrogen storage. It appears that compressed gaseous hydrogen has about 10 times less energy density than gasoline. (Some quick googling found this explanation. It really matters if you're using volume or weight density - by weight hydrogen has a much higher energy density.)

John Tse talked about hydrate clathrates (there's lots of methane at the bottom of the ocean, and it may be a way to store hydrogen ).

One point I found interesting - the potential bewteen a hydrogen molecule and a graphene sheet fits a Lennard-Jones potential quite well (for one particular site or orientation).

Technical note - Doing simulations with DFT is not quite sufficient to get good energies, they used MP2 level theory. (Seems this was repeated in another talk). Question - are MP2 energies good enough, or is doing better just too expensive?

William Goddard presented theoretical evaluations of several approaches to hydrogen storage.

Metallo-Organic Frameworks
Mg nanoparticles
Metal hydrides
Carbon materials (nanotubes, some doped with Li)

In the Mg nanoparticle section, he talked about solving the inversion problem - designed the material to specifications. Their approach was find the potential (I think?) of a desirable atom, and find reals atoms to match.
(There was a similar idea presented by Alex Zunger, trying to find materials that matched specific band-related properties)

Monte Carlo methods in CISE

2006-03-08T22:18:00.000-06:00

Monte Carlo methods are the theme for the March issue of Computing in Science and Engineering. There's a broad overview article, an article covering the geometric cluster method (contains descriptions of the lattice-based cluster algorithms as well), an article examining Markov chains arising from computer science, and an article on sequential importance sampling.

March meeting

2006-03-04T23:01:00.000-06:00

I'll be at the APS March Meeting.

I'm giving a talk in session V27 (on Thursday).

If anyone else going to the meeting wants to meet and talk about QMC, MC, computational physics, or other topics, feel free to email me.

Unbiased exponential estimators

2006-01-18T18:47:00.000-06:00

The penalty method allows a Metropolis random walk to be correct even with the presence of noise in the probability function being sampled.

It eventually dawned on me that this same method can be used to correct any exponential estimator.

Imagine we wish to compute the estimator, exp(x), where x is some sample point that is contaminated with noise. Suppose the variance of the noise is sigma, and we know the noise is gaussian. Then the estimator exp(x-sigma^2/2) will be unbiased.

Free Energy Links

2006-01-15T23:01:00.000-06:00

Radford Neal has written a paper on computing free energy differences. The paper presents the Linked Importance Sampling method, but describes other methods as well. The paper contains references to both the physics literature and the statistics literature, and uses physics language.

If you're really interested in free energy differences, check out Arjun Acharya's thesis, Free Energy differences: Representations, estimators, and sampling strategies (also on arxiv). The review chapter covers several methods, but most of the thesis deals with the Phase Mapping technique, which increases the overlap between phases. It does not eliminate the problem entirely, however, and many of the other methods can be used to deal with the remaining issues.

[Edit 2/27/06: updated links from comments, and revised the entry]

Linear scaling and noisy forces

2005-11-21T20:13:00.000-06:00

This paper by Krajewski and Parrinello,
Linear scaling electronic structure calculations and accurate sampling with noisy forces, has two intriguing ideas.

They compute electronic structure using a new linear scaling method that involves sampling, and thus results in noisy forces. Those forces are used to drive a molecular dynamics simulation of the nuclei (note that one of the authors is Parrinello).

For the electronic structure part of the problem, they replace a determinant in the independent particle grand potential with an integral over some fields. These fields can be sampled using a Langevin equation. This Langevin equation apparently involves only matrix-vector multiplication with a sparse matrix, which delivers the scaling.

The second part (starting in the paragraph before equation 11) is how to deal with the noisy forces in molecular dynamics. Their solution for this is to again use a Langevin equation, and to choose the size of the damping terms to cancel out the noise on average. In this way they obtain correct sampling of the Boltzmann distribution.

Polarization

2005-11-07T22:34:00.000-06:00

Calculating observable properties of systems has been a weak point of QMC methods. The paper, Dielectric Response of Periodic Systems from Quantum Monte Carlo Calculations by Umari, Williamson, Galli, and Marzari remedies some of that deficiency by demonstrating a calculation of polarization using DMC.

"Programming in Math" updated

2005-09-22T00:31:00.000-05:00

The paper on Programming in Mathematical Notation has been updated. (see previous blog entry about Programming in Math)

I've added a page about integration by random sampling. This example adds Greek characters, and a foreign function interface for the random number generator. Also, the translation of several separate sums (defining average, variance, etc) into C++ requires the sums to be merged into a single loop.

[Edit 2/24/2007 - update link]

Cluster algorithm

2005-08-10T18:37:00.000-05:00

In the paper, Cluster algorithm for pairwise-interacting particles, Stephen Whitelam and Phillip Geissler describe an algorithm for creating a cluster of particles to move as a group.

The cluster starts with a randomly chosen particle and grows by adding particles with some probability (based on the Boltzmann factor). Particles that strongly interact with cluster are more likely to be part of the cluster. Then a normal trial move is made of the cluster as a whole.

Importantly, the potential and temperature chosen for this operation are not the same as the potential and temperature of the system. In fact, the cluster size can be controlled by adjusting the temperature (at infinite temperature each cluster is only a single particle). This cluster temperature is varied during the simulation run to sample moves with different cluster sizes.

See the paper for more details about controlling the range of energies involved in the cluster (and the cluster interface energies).

Programming in Math

2005-07-20T18:02:00.000-05:00

In a previous post on Solving science problems, I commented on the process of scientific computation. One of the steps was converting the solvable equations into a computer program.

As part of making that step more transparent, I've been working on method for Programming in Mathematical Notation (works best with MathML capable browser).

The program takes a content MathML file as input and outputs a presentation MathML file for nicely rendered display in the browser, and a C++ file to compile and run.

[Edit 2/24/2007 - update link to example]

Fast variance optimization

2005-05-29T23:03:00.000-05:00

VMC variance optimization usually uses the reweighted variance.

What is reweighting? Consider a set of particle positions obtained from sampling a wavefunction at a given set of parameters. The energy (and the variance of the energy) explicitly depends on the wavefunction parameters by the local energy formula, and implicitly through the particle positions determined by sampling the wavefunction. When optimizing, one would like to vary the parameters without going through all the trouble and expense of sampling new positions. Computing the energy (or variance) change due to the local energy piece is easy - just put the new parameters in and recalculate. Reweighting is the process of correcting the implicit dependence of the energy (or variance) on the wavefunction parameters.

N. D. Drummond and R. J. Needs describe a VMC optimization method in their paper, A variance-minimization scheme for optimizing Jastrow factors, that uses only the explicit dependence of the variance on the parameters (ie, they use the un-reweighted variance). They restrict their scheme to linear parameters in the Jastrow factor, and this makes the resulting parameter surface a quartic. Also, this scheme only requires a single sum over electron positions (the requisite factors can be accumulated during a single run), making the optimization step very fast.

Even though this scheme uses a less accurate approximation for the variance, it appears to give very good results. In fact, the resulting energies for a neon atom are better using the unreweighted variance for optimization than from using the more traditional reweighted variance.

Solving science problems

2005-04-16T07:30:00.000-05:00

What is the ideal programming language for MC simulations?

What if we could design a language expressly for MC and QMC simulations - what would it look like? As a first step, let's put programming in the context of solving the whole problem.

Steps to solve a computational science problem

Define the problem (ask the question)
Write down the equations modeling the system under investigation
Using mathematical transformations and approximations, reduce the equations to a solvable form (numerical analysis)
Convert into a computer program and run it
Undertake verification (are steps 2,3, and 4 correct) and validation (are step 1 and the approximations in step 2 correct)

Developing techniques for managing each of the steps would be a great help. Current programming languages only deal with step 4. Computer algebra systems (Mathematica, Maple, Maxima, etc) can help with steps 2 and 3, but don't really help in making the jump to a large scale program.

Currently the transformations from steps 2 through 4 are done by hand. I know I make many mistakes in this process. It would be a boon to my productivity if the computer could verify that all the minus signs and factors of 2 are correct.

The verification step often involves understanding the effects of various approximations in step 3. If the program is described as a series of transformations on the equations, hopefully changing the approximation would be a simple change to the equations, and the rest of the process would automatically generate a correct program. This would make it easier to test approximations.

This doesn't say much about the ideal programming language, but it does give some idea of the problem to be solved. I'll have more to say about the ideal programming language later.

Programming Language Notes

2005-04-13T19:13:00.000-05:00

On the Sun website, there's an interview with Guy Steele. He talks a little bit about Fortress, which tries to make programs looks more the original problem domain by allowing more mathematical notation. I will be interested to see the results.

Guy Steele is also the principal investigator for the Programming Language Research group. On that page there is a paper titled Object-Oriented Units of Measurement. It's a look the issues involved in making units work in a Java-like OO language.

The last note is an ACM Queue article by Gregory Wilson about Extensible Programming. The article points out that there are advantages to storing the program text in a more structured form, including being able to include diagrams or math notation in the source code.

Whack A Mole

2005-03-24T06:47:00.000-06:00

In QMC, a central role is given to the local energy, defined as (1/psi)H psi, where psi is the trial wave function. If psi is an eigenstate, then H psi = E psi by definition, and the local energy becomes a constant. Just as lower energy trial functions are deemed better than higher energy ones, those with fewer fluctutations are better than those with more fluctutations.

One could imagine looking at the fluctuations in the local energy, and using those to determine where best to adjust the trial wave function to reduce the fluctutations (side note: one would want to use psi^2 times the local energy, since fluctuations are weighted by psi^2 in the integrals). This seems similar to the game of Whack-A-Mole, where one attempts to hit mechanical moles that pop up randomly for a brief period of time. Maybe the method could be called Whack-A-Local-Energy (WHALE).

Unfortunately for my naming scheme (and fortunately for everyone else), the method described already exists, and it's called the Energy Fluctutation Potential (EFP) method. It was first described by Stephen Fahy in "Quantum Monte Carlo Methods in Physics and Chemistry" [NATO ASI Ser. C 525 101, 1999; M. P. Nightingale and C. J. Umrigar, eds]. Further papers describing the method:

Optimal orbitals from energy fluctuations in correlated wave functions, by Fahy and Filippi
Optimization of inhomogeneous electron correlation factors in periodic solids by Prendergast, Bevan, and Fahy
Optimized Jastrow-Slater wave functions for ground and excited states: Application to the lowest states of ethene by Schautz and Filippi

The EFP method involves an iteration that starts with a one-body solver (HF or LDA). The resulting orbitals are used in a VMC calculation, and the fluctuations in the local energy form a extra potential that is fed back into the one-body solver to get new orbitals, and the iteration repeats until a self-consistent solution is reached.

(Name note: this paper by Umrigar and Filippi calls it the Effective Fluctuation Potential method. This paper also gives another reference I didn't put on the above list)

Store, manage, and share your reading list

2005-03-20T16:25:00.000-06:00

The cleverly-named del.icio.us site describes itself as "social bookmarks." The "bookmarks" part allows users to keep track of links online, include a short description, and add an arbitrary number of one-word tags to help categorize the link. With the
"social" part, people's bookmarks are available to everone else.

Here's my list of bookmarks. (I just started, it's not very large)

For the academic world, Richard Cameron created CiteULike. It allows assigning category tags to papers, and has the social aspect of sharing links to papers. Perhaps most importantly, it can automatically import bibliographic information from supported publishers (including arXiv) and export the information to other formats (BibTex, EndNote).

Here's my list of articles. (Once again, very small so far)

The link to CiteULike was found via Corante

I don't know whether either of these will be a useful tool for research, but it seems worthwhile to try them out. When I started grad school in 1993, I went to the library to photocopy articles. By the time I finished in 2000, I went to the library website to download articles. That seems like a large shift, but I suspect it's not just because it makes it faster and more efficient to get the articles - it's the linking, management, and annotations represented by these sorts of applications that will make a bigger change in the future.

MathMath

2005-03-02T21:26:00.000-06:00

A wiki seems like a good way to organize information. Blogs work nice for recording thoughts and links in chronological order, but aren't so good for organizing ideas topically or for revising items. Proto-papers or latex files with background information have been used. The easily edited environment of a wiki seems like it would function well, even if were only used personally (part of the wiki idea is open collaborative editing). The wiki text markup is much lighter than HTML, so it would be easier to deal with than maintaining a set of web pages.

The real issue is integrating mathematics. Mozilla and derivatives can handle MathML (and IE can with a plugin), so that seems like a good route (inserting pictures is marginally okay for converting papers to HTML, but is not as good for a more interactive environment). The only problem with MathML is the rapid onset of RSI if you actually write equations with it.

There are a couple of solutions, both convert TeX or vaguely TeX-like input to MathML. They are itex2mml and AsciiMathML. The second solution seems rather nice, as it uses javascript to convert the TeX-like notation to MathML when the page loads into the browser.

For a wiki, I've been looking into MoinMoin. At first I tried hacking itex2mml into it, and it took quite a few changes to make everything proper XHTML. Then I found the page that describes integrating AsciiMathML into MoinMoin, and that was much smoother.

AsciiMathML also has a nice way of adding new symbols, that doesn't require changing the AsciiMathML.js file.
In keeping with title of this weblog, the Unicode number for hbar is 0x210F. (although it usually only appears once in a paper and is promptly set to "1") The code for adding it is


<script type="text/javascript">
AMsymbols = AMsymbols.concat([
{input:"hbar", tag:"mo", output:"\u210F", tex:"hbar"}
]);
</script>

Chemistry Textbooks

2005-02-20T19:11:00.000-06:00

Jack Simons has a couple of introductory chemistry textbooks on his website. They are An Introduction to Theoretical Chemistry and Quantum Mechanics in Chemistry (with Jeff Nichols).

Stochastic Matrix Inversion

2005-02-10T21:24:00.000-06:00

I ran across this paper, A Monte Carlo algorithm for efficient large matrices inversion. There are a couple of techniques presented, one which is similar to the iterative Gauss-Seidel algorithm.

QMC methods often use a matrix to enforce antisymmetry conditions on a fermionic wavefunction. The determinant is needed for the wavefunction itself, and the inverse matrix is used for the local energy.

I wonder if a stochastic matrix algorithm would be useful in the QMC calculations. My guess is that it would not be (and if there is an efficiency gain, it might be offset by the added complexity), but it would be interesting to try, if for no other reason than ideological purity (if you're going to use a stochastic method, might as well make it as stochastic as possible!)

Other interesting questions arise. The locations of the all-important wavefunction nodes will now be fuzzy. What consequences does this have for a fixed-node algorithm?

Inverse matrices also appear in Sorella's Stochastic Reconfiguration method for optimization. Another use?

On a side note, I tried looking for background on Google, but searching for "stochastic matrix inversion" returned results relevant to inverting a stochastic matrix rather than stochastic methods for inverting a matrix.

Nifty math and physics applets

2005-01-12T20:22:00.000-06:00

Paul Falstad has a page of math and physics applets. (found via Alan Aspuru-Guzik's weblog)

Get your units correct

2004-12-15T06:42:00.000-06:00

Check out Frink. It's a language designed to ease calculations involving units and unit conversions. I haven't tried it, but it looks interesting. The examples are amusing, and the data file for all the conversions (units.txt) makes interesting reading. Also interesting is Alan Eliasen's presentation at the Lightweight Languages 4 conference.

Hierarchical VMC (was Probabilistic Functions III)

2004-11-21T22:30:00.000-06:00

A couple comments from the last post.

I had wondered whether it was possible for VMC to get as good of results as DMC. In theory, of course, the answer is "yes". Provided there is enough variational freedom. This is the hard part.

On to the second point.

I had proposed a variational wavefunction consisting of the sum of gaussians. The centers and widths were then the parameters to be optimized. The starting center points were to be generated by sampling from a more traditional wavefunction. Why bother? Why not start from a uniform distribution?

My speculative answer is that it would be horribly inefficient, and very likely to follow wrong paths (since the wavefunction is so "floppy"). So starting from a traditional wavefunction is a form of accelerating the convergence of the optimization procedure.

The idea seems to be one of hierarchical VMC, or of a hierarchical wavefunction that can be optimized in steps. The coarse parts are optimized first, and then the finer details are added later. Vaguely like a wavelet decomposition or a multigrid method deals with spatial resolution.

It still remains to be seen if the sum-of-gaussians wavefunction survives contact with reality.

Probabilistic Functions II

2004-11-18T22:11:00.000-06:00

Last time I looked at the sequence of values obtained from sampling a function. I was trying to understand what would happen if used a sum of gaussians, each centered at one of these values.

I was confused because I had a set of sample points, but also a function using those points that could be considered a variational trial function in its own right.
The width of gaussians can be taken as a variational parameter, and optimized.

The idea for a MC scheme would be to generate a set of sample points from a trial function. Then the diffusion algorithm would be stepped forward in time to get a better approximation to the wave function.

But we can go further and consider all the center points as variational parameters as well. And we don't actually care about the diffusion part - it is only a device to relax the wavefunction to the ground state. We could use any optimization method instead. Now this form for the wavefunction is likely to have lots linear dependencies. Some optimization methods may have problems with this.

One could simply do a random walk - try moving a point, evaluate the energy (with VMC), and accept or reject the move based on the energy difference. Could one get DMC-like quality from VMC?

I can see a couple possible problems, so far. Getting the boundary conditions and cusp conditions correct may be difficult. The "basis" functions don't need to be a set of gaussians - just positive functions with a center point and a width. Also, there may be so much variational freedom that convergence takes a long time.

Probabilistic Functions

2004-11-17T22:34:00.000-06:00

A function, f(x). We can treat it as a black box - feed it x, and get back f(x). But what if we have different sort of box - one with a big green button. Push the button and it spits out an 'x' value. The only way to ascertain f(x) is to keep track of the density of the values of 'x' that come out of the box - where f(x) is large, there will be lots of x's, and not so many where f(x) is small.

If we want to compute a normalized average over a weight function f(x), this output is exactly what we want. Now suppose we want to compute derivatives of f(x) (the energy in QMC calculations). In VMC, we can open the box and there is a analytic form for f(x) sitting there that we can take derivatives. But in DMC, the situation may not be so simple - f(x) may be determined solely by the distribution of x's. (I was trying to figure out how to compute the energy when doing DMC w/o a trial function)

One way to compute the derivative of f(x) is to keep a histogram of x's and take the derivative after the fact. However, this solution seems noise prone, and not likely to scale well in multiple dimensions.

Each sample is a delta function. Integrals over the function become integrals over a sum of delta functions, centered around the sample points. What if we replaced the delta functions by a gaussian? (since that is delta function in the limit of vanishing width). Each point then represents a probability distribution with a finite width.

The first attempt at representing the ground state wavefunction of an infinite square well failed miserably - I could get any value of the energy I wanted by varying the width of the gaussians. The problem is the boundary conditions - the wavefunction must vanish at the edges. So then I tried using image gaussians (outside of the square well) to force the appropriate BC's. That worked much better. So well, in fact, that the value of the energy is quite good even when the center points are drawn from a uniform distribution, although I suspect this is largely a feature of simplicity of the potential).

So what now?

With this sort of representation, I'm wondering if it's possible to use a forward Euler scheme to propagate the distribution forward in time (ie use a first order approximation to the time derivative in the diffusion eqn.). Either symbolically (ie, complicated formulas related to the original points) or by sampling (ie, generate a new set of points).

Normally, one thinks of using orthogonal function expansions, since linearly dependent functions don't add anything new to the function. This sum of gaussians is definitely not orthogonal (well, they become more so at small widths), especially since the points are randomly distributed according to some distribution. Hmm. This contrast seems relevant somehow. The gaussians are all positive, and orthogonal functions usually have negative regions - I don't know if this is relevant or not.

Comments for Everyone. A New Blog.

2004-11-11T20:58:00.000-06:00

I turned on comments for everyone, not just registered users.

Also, feeling the need to speak of things other than QMC, I've started a different blog for such topics. (I will continue posting here, I just wanted a place for non-QMC or non-MC related items).

[Edit - removed link to personal blog. Email me if you really want to know]

Atomic and Molecular orbital applets

2004-11-11T20:17:00.000-06:00

Here's a course page that contains some atomic and molecular orbital applets. Look under the "Dry Laboratory Experiments", "D3" for atomic orbitals and "D4" for molecular orbitals.

DMC timestep and equilbration

2004-11-04T20:44:00.000-06:00

I've been reading Chebyshev and Fourier Spectral Methods by J. P. Boyd.

In Chapter 12, p. 227, he points out that if we are only interested in the steady state, errors in the decay rate don't matter.

Chapter 13 talks about operator splitting, which looks very similar to how Diffusion Monte Carlo gets the short time approximation to the propagator. Section 13.4 talks about consistency - basically timestep error in the steady state.

I had the realization that one could use large timesteps in DMC for the equilibration part, and smaller timesteps once equilibrium was reached. More drastic population control measures are probably needed as well if large timesteps are used (at least I always had trouble maintaining a stable population of walkers if the timestep is too large).

This particular idea isn't terribly profound, but I'm finding it very enlighting to read about other methods for solving differential equations. I also like the presentation level of this book - lots of practical, concrete considerations as well as descriptions of the consequences of various equations.

Qumax - QMC code

2004-11-03T22:37:00.000-06:00

I signed up for Google Alerts a while back, and today I got a notice about Qumax, a Quantum Monte Carlo code.

I downloaded the code, but haven't tried compiling or running it yet.

Fermion Sign Problem Officially Hard

2004-10-27T23:26:00.000-05:00

It now official: the fermion sign problem is hard. Matthias Troyer and Uwe-Jens Wiese present a proof that the fermion sign problem is NP-hard. They do this by mapping the quantum problem to a spin glass model, which is known to be NP-complete.

Do I believe it?

Hmmm. The results presented deal with a discrete problem. I'm not sure if the results apply to continuous systems as well or not. The whole notion of NP-completeness and combinatorial complexity deals with discrete systems. For continuous systems, information-based complexity (IBC) might be a better way of analyzing the problem. Note that the "fruit-fly" of IBC research is multidimensional integration.

However, the distinction between continuous and discrete might not be important. For instance, with global optimization, finding the locations of the local minima is "easy", but then one is left picking between many different discrete local minima to find the global minima (which makes it a potentially "hard" combinatorial problem - assuming the number of local minima grows exponentially with problem size).

I'm still perplexed as the origins of the difficulty of the FSP (actually, I don't understand the origin of difficulties for most NP problems.) So a quick look into the form of the integrals is in order.

We start with wanting the integral I = int f(x) dx, where x is really a vector of high dimension. Monte Carlo fits the bill nicely, having error bounds independent of dimension. But f(x) is strongly peaked - only small volumes contribute to most of the integral. To be efficient, we need importance sampling. The usually trick is to split f(x) into a product: p(x)g(x), where p(x) >= 0. Since p is positive, we can interpret it as a probability and use various sampling methods to get int p(x) g(x) dx / int p(x) dx.

We can get away with this in QMC by assuming that p(x) = psi^2(x), which is positive. But in Diffusion Monte Carlo, p(x) = psi(x) phi(x), where psi(x) is a guess and phi(x) is the exact ground state wavefunction. Now psi and phi will not have zeros at exactly the same location, so p(x) will have some negative regions.

The trick for dealing with negative p(x) is to split it once again into: sgn(p) abs(p) and sample abs(p), and move sgn(p) over with g(x). The integral now looks like
[int abs(p) sgn(p) g(x) / int abs(p)]/[int abs(p) sgn (p) / int abs(p)]

Here's where I get lost. One source of the problem could the denominator - since it is the average of a signed quantity, it could be small (at least relative to the errors) and magnify errors in the numerator? But the paper says that both the numerator and denominator have errors that grow exponentially in system size. So the denominator is an irritant, but not the root cause.

The error of the average sign (again from the paper) is related to the exponential of the free energy difference between the fermionic problem (with p(x)) and bosonic problem (with abs(p(x)). So the fermionic free energy is necessarily smaller because of cancellation induced by minus signs. But I still don't see how this makes the problem hard.

Why not add a sufficiently large constant to p to make it postive? Now there are no minus signs anywhere. But the essential features of the integrand haven't changed - so it can't be a "sign" problem. I guess I view the complexity of integration as coming from the unpredictable fluctations of the integrand - wildy fluctating integrands are "harder" to do than smooth integrands. Adding a constant doesn't change this.

Perhaps the antisymmetry requirements produce a function that has nearly equal positive and negative parts. Each one separately has a small relative error, but the nearly exact cancellation causes great magnification of the error (a well known numerical analysis bugaboo). Thus, the experiment of adding a constant doesn't change anything - it's added equally to both parts, and cancels exactly.

Lennard-Jones Server, Part II

2004-10-27T19:33:00.000-05:00

Some more thoughts, in a Frequently Asked Questions (FAQ) format. (Actually, no one has asked any questions, so maybe it should be Frequently Unasked Questions. Or maybe not, given the acronym. Perhaps I should call them Frequently Anticipated Questions.)

Q. Why remote computation - why not just download a program to the users computer?
A. Have you ever downloaded software and tried to obtain all the dependencies? Hopefully, the client program would be much simpler. Other ideals are bug fixes and enhancements would take place on the server and automatically be picked up by the users (although this would require versioning to meet reproducibility requirements for researchers). The general idea of web services is somewhat of a holy grail for business types - getting paid every time someone uses the service, or getting paid on a subscription basis. Which is why word processing software is unlikely to take off as a web service any time soon. Getting someone to pay for a Lennard-Jones server is also unlikely.

Q. Why not just have a web page? Like Google, Amazon, etc.
A. A web service allows interfacing to these services under program control, which allows your computer to buy stuff on Amazon when it gets bored.

Q. Why a single component Lennard-Jones system? It's not exactly a hotbed of research.
A. No, but it is a testbed. New techniques are often tested on simple systems, and benchmark results would be useful. More importantly, it is a simple first test of the simulations over web services idea. I can see expanding to more complicated classical systems and QMC later.

Q. Can you explain more about accumulating a database of properties?
A. Yes. Aggregating data is a noble and useful goal, but it requires careful vetting of the input data. The idea behind a Lennard-Jones server is that users would only input physical parameters, not raw data. The data aggregation could then be automated, since the data (ie, the output of the program) should be of good quality (or at least of known quality).

The particular example I have in mind is computing free energies and the phase diagram. Computing a grid of pressures or energies at a number of state points is relatively easy, but computing free energies and such is harder - it can require integrating over a number of state points (if you use thermodynamic integration). Hopeful as more state points get put into the system, the free energies get more accurate. Calculations done during the idle time could be planned specifically to lower the errors on edges of the phases.

Now such a database will likely be inferior to targetted research trying to answer a particular question (ie, accurate freezing points, etc), but it seems like a good source of background information or a starting point. Also, the background or baseline information would be in electronic form, and as such be more useful than in a paper. (Not that papers aren't important - the raw data needs to be intrepreted and understood)

Q. Who would be likely users for such a service?
A. On the computation side, companies and researchers that use the results of simulations, but don't care about the details of the simulations. For researchers in the particular field the web service covers, it might not be so useful (the whole "not outsourcing your core competency" thing), although it would be useful for benchmarking and comparison of results.

On the data side, having a good collection of properties for the desired system would be useful for everyone.

Q. What would you do if it were the late 90's?
A. I would run out and register eSimulations.net.com and/or eWavefunctionsRUs.com, write a buzzword compliant business plan about the wonderful opportunies in B2C, B2B, R2B and R2R markets (the last two, made up on the spot, are Reseacher to Business and Researcher to Researcher.), and how actual experiments in the real world are passe, and the new e-world of the internet will make the real world obsolete. And try to get lots of VC money and go public with an overhyped IPO. Ah. For the good old days.

Comments enabled

2004-10-24T22:44:00.000-05:00

I've enabled comments. Just for registered users for now.
And I've started adding titles to the posts. Enjoy

2004-10-23T21:21:00.000-05:00

Lennard-Jones Server

If I could write blog entries as funny as Rory Blyth, I'd be famous. Well, maybe not. I'd have to draw cartoons, too. But to really be famous, I'd have to gratutiously link to Robert Scoble. (And get him to link back to this blog, the really hard part.) (For those who have non-computer-geek lives, Scoble wrote an article on how to get your blog noticed.)

The vaguely relevant part of that last paragraph is Rory's comments about XML Dev Con 2004. Some of the talks related to Web Services.

Web Services - what are they? Zef Hemel has a good description of Service Oriented Architecture, which basically means using web services (ie, Remote Procedure Calls over port 80) to construct an application.

That leads me to the Grid, of which web services seem to be a part. (But don't get me started on the popular analogy between the power grid and the computational grid. It's seductive and appealing, and perhaps even useful at a very abstract level, but the details break down badly if you push it very far (see the Scientific American article, among others.)

The grid concept seems to encompass four different things (as far as I can tell)

Connecting geographically (and/or institutionally) separated supercomputers to work on a single calculation.
Harvesting idle cycles on desktop PC's
Running specific computations on remote machines (this is web services)
Running arbitrary computations on remote machines (like the recent Sun announcement, but other companies have offered similar services previously)

Concerning items 3 and 4, I suspect web services will be more popular then running arbitrary computations remotely. If you have enough computational or data requirements to entertain thoughts of rent-a-FLOP (or rent-a-MIP), you will care about the environment and how well the code is optimized for the target architecture (or at least how predictable the run time is). But caring about such things makes buying "generic" computing power harder, or at least more complex.

But web services are (?) about hiding those nasty details. The provider can optimize the code for the target architecture or other details, and the user need not worry. It may be useful to divide the services into two kinds, one that requires lots of data, and the other that requires lots of compuations. For data-heavy services, the provider has lots of specialized data that is updated frequently (think Google, Amazon, Ebay, or stock market info), and delivering the data + program to the user would impractical. (versus something like TurboTax, which contains lots of specialized knowledge that changes yearly, but is practical to deliver the data + program to the user's computer. Although it can be and is offered over the web as well). For computation heavy services, optimization for particular hardware, etc could be considered specialized knowledge as well (making the divisons less distinct), but I still think the division may be useful.

As a way of understanding web services, and grid concepts, I was thinking of a Lennard-Jones server. A user could request properties of a single component system interacting via a Lennard-Jones potential for a given temperture/density/system size, the server would run the simulation and return the energy,pressure, g(r), and whatever else. (It would require limits on the system size to keep from completely overloading the server).

The server could have other features, such a keeping a database of run results, and returning the simulation results nearest the requrest parameters. Although I'm not sure this would work so well if one extended the idea to more complicated systems with larger parameter space.

In some ways, it's like a database of material properties, except that the properties are all generated by computation. So the database could be improved by further computation - if the server is idle it could work on reducing errors in various quantities. Although this might wreak havoc on reproduceability for users - good versioning would be needed.

2004-10-16T10:51:00.000-05:00

3D Function Plotter

In the continuing quest to visualize wavefunctions, I've put up the first version of a 3D function plotter.

2004-09-29T23:55:00.000-05:00

How to Integrate a Derivative

In the previous post, I had tried to use derivatives to improve the estimate for a single sample, but that turned out not to work very well. That was an approximation. To use the derivative exactly, the definition is f(x) - f(0) = int_0^x f`(y) dy.

Does this have a practical application? After all, if we do this on the example in the previous post, we go from a single integral to a double integral - it looks like we're making the problem harder! Ah, but this is Monte Carlo. The convergence rate is independent of dimensionality. So adding more integrals doesn't necessarily make the problem harder.

The final error does, however, depend on the intrinsic variance of the integral. In this example, this variance gets much worse, so making this transformation is not a good idea. In general, derivatives will vary more than the original function, so this technique is not likely useful in general. Unless the transformed problem has some other desirable features, like better importance sampling or something.

Finally, I should note that something similar to this is used in thermodynamic integration to get free energies.

2004-09-27T18:45:00.000-05:00

Sometimes it doesn't work

The other day I had an idea about trying to improve MC samples. The Monte Carlo method for evaluating an integral is the average of a set of samples. Each sample can be considered a very crude approximation to the integral - each sample approximates the integrand as a constant. What if we could make each sample a better approximation, say by using a the derivative information at the sample point?

Another consideration is that the intrinsic variance of the integral (int [f(x)-ave(f)]^2) controls the accuracy of the MC approximation. If we could smooth the integrand, and lower the variance, the MC integral would be better. Perhaps expanding the integral in a series around the sample point and integrating would help smooth it?

Then reality struck. I tried it on v(x)*exp(-v(x)/T), with v(x) = 1/x^12. I made the expansion and integrated in an interval from [x-h,x+h], and graphed the result - it was clearly *not* smoother than the original function. It was also clear that the expansion was an approximation, and only correct in the limit of small h.

Upon further reflection, I realized a few things. One is that each MC sample may be a crude approximation to the integral, but it has the important property of being unbiased (ie, all the errors will average to zero). Any improvement in the sample value needs to retain this property.

Secondly, smoothing would be most valuable when the integrand is rapidly changing. And this is just the condition where the Taylor expansion is likely to be the least accurate. I suppose that any local smoothing of the integrand will suffer from this sort of problem.

Sometimes, it doesn't work. Ah, but that's what research is all about.

2004-09-04T14:35:00.000-05:00

Visualizing Wavefunctions

I've started learning VTK, Python, and Tk in order to look into visualizing wavefunctions. So far, I've made an application that plots a function of x and y. The user can change the plotted function in an edit box. Get the 2D plotter here.

I've also got a prototype application that plots a function of 3 variables, and lets the user choose a slice (a plane) through the displayed volume, and the application plots the 2D slice in a view that looks very much like the 2D plotter. But it's not quite ready for public consumption yet.

2004-07-29T00:39:00.000-05:00

Welcome to the Metropolis Disco - How do you move your particles?

Out on the dance floor tonight we have the spotlight on two papers on improving the efficiency of Monte Carlo moves. But first, a brief summary of the results of Peskun's theorem (Both papers mention it.). The Metropolis algorithm takes a trial move and accepts or rejects the move with some probability. If the move is rejected, that location is used in the average again. Rejected moves will increase the autocorrelation, and so increase the variance of the resulting average. Peskun's theorem validates this intuition, under certain conditions. So the result is that, just like dancing, moving is better than standing still.

The first paper is Delayed Rejection Variational Monte Carlo by Bressanini, Morosi, and Tarasco. The problem is that for best efficiency, moves near an atomic nucleus should have small step, and moves far away from a nucleus should have a large step. If a move is rejected, it is likely that the step size was too large. One could imagine a dynamic or adaptive step size method that lowers the step size when there are lots of rejections and increases it when there are lots of acceptances (*). This paper does something slightly different - when a move is rejected, it is retried with a smaller step, but the original rejected move is not made part of the average. The acceptance probability is adjusted to take this into account. So there is an increase in the amount of time it takes for a move, but that is made up for by a larger decrease in the autocorrelation time.

The second paper is Optimal Monte Carlo Updating by Pollet, Rombouts, Van Houcke, and Heyde. They work with discrete systems, and introduce a a procedure called "Metropolizing" - which is to take a move rule (eg, heat bath) and trasform it to a new, more efficient transition matrix. This process can be applied repeated until there is only one non-zero diagonal element, resulting in a very efficient transition matrix (For single site updating rules).

(*) I'm not entirely sure a simple dynamic step size method would satisfy detailed balance. I think it would, but that should be checked.

2004-07-21T15:46:00.000-05:00

More Graphic Card Computations

There's a site, GPGPU (General Purpose computation on Graphics Processing Units), dedicated to interesting things one can do with a graphics card. (Well, there's always graphics, but using things for their intended purpose can get boring :)

Last fall, I wondered about the possiblity of doing Monte Carlo on the graphics card.

On the GPGPU site there's a link to a paper about doing Ising model and percolation simulations on the graphics card (Benchmarking and Implementation of Probability-Based Simulations on Programmable Graphics Cards). The work is primarily for proof of concept and benchmarking, with an eye towards lattice QCD simulations eventually. Hmm. Play a killer game of Quake *and* uncover fundamental secrets of the universe - what a combination!

I can imagine students requesting the latest "Ising co-processors" so they can "make progress on their thesis" when Doom 3 comes out. (Err, maybe they need to leave off that last part when making the request).

2004-05-12T23:47:00.000-05:00

PIMC live

Check out the Path Integral simulation applet on John Shumway's page. You can change the temperature and see the paths localize at high temperature and spread out at low temperature.

2004-04-16T01:03:00.000-05:00

Pseudorandom number generators and the stock market

As correlations are discovered in psuedorandom number generators used in simulations, the correlations are removed by creating better generators. So as time marches on, the state of the art PRNG's should be getting more and more random.

As correlations are discovered in the stock market, they are removed as the resevoir of correlation is drained for its profit. So the stock market should be getting more and more random as time goes on (especially as faster computers and more sophisticated statistical techniques are used)

2004-03-11T22:59:00.000-06:00

Experimentally switching between liquid and solid phases?

Looking at Grochola's free energy path, I thought it was a nice computational technique.

The March 2004 issue of Physics Today has an article on controlling atoms and ions in optical traps (it's actually about using these for quantum information processing). On page 41, it describes putting the atoms in a optical lattice. The depth of the wells and the hopping matrix elements can be controlled by varying the intensity of the lasers creating the lattice.

This seems like one could experimentally switch between solid and liquid states without going through a phase transition.

There are some differences - I don't know if they are important or not. First, the experimental setup is a quantum system, Grochola's paper describes a classical system. Second, the confining potential is different. In the experimental setup, atoms are confined to any lattice site, but are not tethered to a particular lattice site.

The Physics Today article mentions a Mott insulator transition when the confining potential becomes strong enough. This transition seems due to the superfluid (quantum) nature of the liquid state. If the experiment were performed on a non-superfluid system, this transition would not happen ??

Update (3/16/2004): I forgot the final stage in the path - turning off the confining potential to get a stable solid state. I'm guessing the lattice spacing of the optical lattice is at least the laser wavelength or greater, so it's much larger than the lattice spacing of a normal solid. But some atoms can have their scattering length tuned - could it be made large enough to be comparable to the optical lattice spacing? In that case the atoms would behave as hard spheres (most importantly, there would likely not be enough attraction to keep the solid together), so an external potential would be required. Most of these setups seem to generate a harmonic confining potential (a box-like potential seems much harder to generate??), so the question is - what is the behavior of hard spheres in a harmonic potential?

For this system, there are two parameters - the number of particles and the strength of the external potential. I suspect there would be three regimes - one all liquid, one with a solid center surrounded by a liquid exterior, and one would be all solid. To take the thermodynamic (large N) limit, it seems that one would need to increase the particle number and decrease the external potential strength simultaneously to keep the central density or pressure constant (otherwise the pressure/density at the center would increase with added particles). Or maybe such a limit doesn't exist?

The Way of Antisymmetry

2004-02-11T23:59:00.000-06:00

The traditional QMC wavefunction is the product of a determinant of single particle orbitals and a Jastrow factor (which contains the electron-electron correlation). The determinant ensures the required spatial antisymmetry of the wavefunction, and is computable in O(N³) time (where N is the number of electrons). But it is not the most general form possible.

The most general form is obtained by explicitly antisymmetrizing: summing the values of a function with all possible permutations of the electron coordinates (and multiplying by the sign appropriate for the permutation). For example, with 2 particles and a function f(r₁,r₂), the antisymmetrized form is f(r₁,r₂) - f(r₂,r₁).

Unfortunately, the number of terms is N!. (The determinant form is obtained by applying this procedure when the function has the special form of the product of functions that each depend only on a single electron coordinate - these functions are the single particle orbitals. In the two particle example, f = f₁(r₁)f₂(r₂). The antisymmetrized form is f₁(r₁)f₂(r₂) - f₁(r₂)f₂(r₁), which is the determinant of a matrix.)

Dario Bressanini, Massimo Mella, Gabriele Morosi, and Luca Bertini have several papers using this form of the wavefunction for small molecules, and they get very good results. The papers can be found on Massimo Mella's publication list (edit 12/3/2005 - removed nonworking link). They are (this may not be a complete list):

Many-electron correlated exponential wavefunctions. A Quantum Monte Carlo application to H2 and He2+

Chemical Physics Letters 240, 566 (1995)

Paper available at this Citeseer page
Nonadiabatic wavefunctions as linear expansions of correlated exponentials. A Quantum Monte Carlo application to H2+ and Ps2

Chem. Phys. Lett. 272, 370 (1997)

Paper available at this Citeseer page
Linear expansions of correlated functions: a variational Monte Carlo case study

Int. J. Quantum. Chem. 74, 23 (1999)

Paper available at this Citeseer page
Explicitly correlated trial wave functions in quantum Monte Carlo calculations of excited state of Be and Be-

J. Phys. B, (2001)

Not all N! terms will give a significant contribution to the result. Is it possible to use techniques similar to O(N) approaches to reduce the number of terms? Is it possible to sample terms from the sum? These might have trouble near the nodes, since the wavefunction vanishes due to the cancelation of all the terms.

2004-02-11T14:39:00.000-06:00

DMC of porphyrin

An all electron calculation of a large molecule (large by QMC standards, anyway):
Quantum Monte Carlo for electronic excitations of free-base porphyrin by Alan Aspuru-Guzik, Ouafae El Akramine, Jeffrey Grossman, and William Lester. The porphyrin they used has the formula C₂₀ N₄ H₁₄, so 168 electrons.

2004-02-04T16:10:00.000-06:00

Solid vs. Liquid, The Free Energy Battle

Determining the free energy difference between the solid and liquid states is essential for accurately mapping the phase diagram of a substance. Typically, one determines the free energies of the solid and liquid states through separate thermodynamic integrations, where each integration has one end point connected to a state with an analytically known free energy. In the solid state, the Einstein crystal (springs attaching particles to lattice sites) is often used as a starting point. In the liquid state, the low density gas can be used as an end point. Other more complicated scalings of the potential are also possible for the integration path.

A direct integration path between the solid and liquid states would be useful.
Gregory Grochola has published one such path in "Constrained fluid lambda-integration: Constructing a reversible thermodynamic path between the solid and liquid state"(JCP 120, 2122).

The basic path is to start with the liquid state and turn off the intermolecular potential. Then turn on a potential constraining particles to lattice sites (He used a Gaussian potential). The last step is to turn off the lattice potential while turning the intermolecular potential back on.

He mentions at the end that this method is not very computationally efficient. I would expect using Bennett's method for free energy differences would cut down on the number of intermediate states needed.

Questions

Compare and contrast this path with the phase switching Monte Carlo technique (there's a reference in the article)
Thermodynamic integration requires there be no phase transitions on the path. Is it immediately obvious that this path has no transitions? Under what conditions, if any, would a phase transition occur on the path?

2004-01-24T16:11:00.000-06:00

Hard Spheres on the Move

Cluster algorithms increase efficiency for systems on a lattice by flipping large numbers of connected spins at once. In "Cluster Monte Carlo algorithms", Werner Krauth summarizes work on how such algorithms can be used in a continuum hard sphere simulation.

The basic idea is to rotate particles by &pi around a randomly chosen axis. (The key observation is that this transformation is its own inverse). A practical algorithm (the pocket algorithm) chooses a collection of particles that can be rotated without any hard spheres overlapping. It starts with a single particle and rotates it. If it overlaps with any other particles, they are added to the collection. Then the new particles are tested in the same way, and the collection is built up of all the particles that need to be rotated to prevent any overlaps.

At high densities, this algorithm will fail because the pocket cluster will contain all (or most) of the particles in the system, and simply rotate the entire configuration. Despite this limitation, there are systems where these moves are useful (see the paper for some examples).

On another subject, suppose you wanted to keep the trial box size fixed, and see how the dynamics varied with density (which would be useful if one were studying the actual MC dynamics) At higher densities, the efficiency would get quite bad. Hirosi Watanabe, et al tackle this problem in A rejection-free Monte Carlo method for the hard-disk system. They use geometry rather than random sampling to determine the allowable move area (and acceptance probability). To actually choose a point, they use a rejection method (it's not clear why this method is called "rejection-free", then. Maybe because the rejections are not essential for maintaining detailed balance, and rejected moves don't need to be included in the average.). At the step of choosing a new point, the area they choose from is greatly reduced from the original trial area.

2003-12-12T00:46:00.000-06:00

Reducing the Variance

Assaraf and Caffarel have published a couple of papers (JCP 119, 10536 and JCP 113, 4028) about reducing the variance for general estimators in QMC. The VMC energy has the well known zero variance property - the variance goes to zero as the trial wavefunction approaches an eigenstate. Through a "renormalization" , the authors show how to extend this property to quantities other than the energy.

The basic approach seems to be to add a term (with zero average) to the bare estimator. Then that added term can be optimized to reduce the variance. Their derivation starts with the Hellmann-Feynman theorem, but this seems confusing to me.

Why not start with the new estimator, O tilde = O + B. ("O tilde" and "O" are notations from the paper. "B" is my own notation for representing all the added terms) Then we demand < B > = 0, and use a Lagrange multiplier to enforce it while optimizing the variance (<(O+B)^2> - < O+B>^2)?

Their paper gives a particular form for B, that ensures the < B>=0 constraint is met (or at least that < B > vanishes as the trial wavefunction approaches the ground state).

Can this approach be used to reduce the variance in classical Monte Carlo as well? In principle, B could optimized in that case as well - the hard part is finding a special form as in Quantum Monte Carlo.

Added thoughts (12/14/2003): The primary application of this method is to forces, which are plagued by an infinite variance in QMC. This looks a like good way to improve computations of the force.

The paper looks at average values and variances as functions (or at least the polynomial order) of delta psi (the diference between the exact ground state and the trial wavefunction). For general (non-energy) estimators, the bias (or error) in the average value is linear in delta psi, and the variance is constant. Using their procedure, the error is quadratic and the variance is linear. It appears the added terms (B) cancel out the leading error terms in the estimator. The question is then, why stop at the leading order? Can a similar process by applied to cancel out the quadratic term in the average ??

2003-12-04T18:17:00.000-06:00

Free Energy Chronicles

Free energy calculations are annoyingly difficult. Computing averages with Metropolis is relatively easy,
but the free energy is hard to get. The usual method is thermodynamic integration over a path from a state of known free energy to the desired state. (Maybe it's the ensemble averages that are anomalously easy, making other quantities seem difficult by comparison) In addition, some other aspects about free energy calculations are bothering me as well.

The typical example for explaining Monte Carlo integration is to compute the area of an irregularly shaped field or pond by throwing rocks at it. But then computing the area (volume, free energy, whatever) is actually quite difficult in practice. Of course in this example, one has to compare the rocks falling in desired area to those falling in an enclosing area that can be computed analytically, so it's still computing the difference between two states.

Another issue that bothers my intuition is that if I were to wander around in a field I would have some idea how big the field is based on my wanderings. In other words, why can't we get the area from properties of the random walk. This could simply be pushing the analogy too far. A person wandering is not undergoing a truly random walk. They would likely walk in single direction until hitting an edge and then change direction. So imagine riding an out-of-control Segway around a field - would you still get a sense of the size of the field.

Additionally, a person is able to see the edges. So imagine being in a corn field (I am in central Illinois, after all - others can imagine a forest, with very soft trees) on an out-of-control Segway. Would it still be possible to get an idea of how big the field is based on the path taken and how often you bump into the edge?

So I'd like to look into some of the properties of random walks - how much information can be extracted. Are there theoretical limits? Can we change the random walk or come up with different estimators that would let us compute the area?

In order to investigate concrete but simple systems, I'm thinking of a series of systems to make the problem increasingly realistic for physical systems. The simplest and easiest to think about the area of an irregular region. More realistic, the area of the region is much less that the area of the enclosing square, and the region is long and snaky. Then go to higher dimensions, and eventually to a system of hard spheres. To be more generally applicable, the next step would be a general potential (think of a hilly field and we want to know the volume of soil is in the field)

2003-11-10T20:54:00.000-06:00

Interesting notes

While searching for Monte Carlo information on the web, I ran across Jonathon Goodman's Monte-Carlo Course Materials. One interesting part was Alan Sokal's introductory notes near the bottom (see also this CiteSeer link). I liked these notes because they covered a practical method for computing the autocorrelation time.

2003-11-02T23:43:00.000-06:00

Optimization

Recently there were two papers in the Journal of Chemical Physics related to QMC optimization. The first, Geometry optimization in quantum Monte Carlo with solution mapping: Application to formaldehyde, is concerned with optimization of nuclear positions using only energy information, and not gradient (force) information, since it is noisier and harder to obtain. They use statistical experiment design techniques to pick coordinates at which to perform the energy computations, and fit that data to a quadratic form. The the geometry was then optimized using that fitted form.

The second paper is concerned with optimizing the one-body part of the wavefunction: Transcorrelated method for electronic systems coupled with variational Monte Carlo calculation. The authors use a transcorrelated Hamiltonian (some of the correlation has been moved into the Hamilitonian), and compute the SCF solution for the one-body orbitals. They use regular VMC to optimize the Jastrow factors. Note that the two calcuations are coupled, and to reach a solution the procedure is to alternately optimize the one-body and Jastrow factors.

2003-10-08T15:09:00.000-05:00

Monte Carlo on the Graphics Card

A recent short paper on The Graphics Card as Streaming Computer (see related presentation slides) got me wondering if Monte Carlo could be programmed onto a graphics card.

The proposed Protein Explorer from Japan, using a GRAPE type architecture, follows a similar programming model.

One limitation is the precision of the calculations. Single precision at best, some cards only have 8-16 bits. The penalty method or related methods may be useful for correcting errors due to finite precision. (See also Ball, Fink, and Bowler)

Another limitation is the size and complexity of the code that can execute on the graphics card. This might be more useful for simple classical Monte Carlo simluations that QMC.

I suspect this falls into the category of "cool ideas that would be fun to play with, but spending the time researching better algorithms is probably more effective".

2003-10-02T00:45:00.000-05:00

Of Wavefunctions and Wikis

and mailing lists (but that didn't alliterate nicely in the title).

There's a new QMC Mailing list in town. Be the first on your block to sign up.

Also there's a QMC Wiki (A wiki is an easily editable set of community web pages) (No link yet - you'll have to sign up for the mailing list and look at the Sept. archives to get the address)

Both of these are brought to you by Alan Aspuru-Guzik.

2003-09-11T22:53:00.000-05:00

Linear charged particle simulations

Today's paper is by Joerg Rottler and A.C. Maggs, entitled "A Continuum,O(N) Monte-Carlo algorithm for charged particles". There are some earlier papers in the references where this method was developed for lattice systems.

As near as I can tell, the algorithm works by writing the total energy in terms of the electric field rather than the electric potential. The electric field is put on a grid. Particle moves then only require local updates to the electric field grid.

The electrostatic interaction is generated as a result of minimizing the energy functional (functional of the electric field). The algorithm they give is supposed to recover the proper interaction by a Monte Carlo means. This is the part I don't quite understand yet.

It seems like a good example of exploiting Monte Carlo to reduce the complexity of a problem by solving it on average.

2003-07-23T23:31:00.000-05:00

Linear scaling QMC

"Linear scaling of the local energy in quantum Monte Carlo" JCP 119, 1307 by Manten and Luchow.

Cutting off the orbitals and using sparse matrix computations for the determinant, and cutting of the Jastrow factor are the time consuming parts.

One interesting point - they use a topological cutoff for the orbitals rather than a strictly spatial cutoff.

This is the scaling of the local energy with the number of electrons - I'm curious how the correlation time varies with system size.

2003-06-18T18:06:00.000-05:00

50th Anniversary of the Metropolis Algorithm

Last week I attended a conference on the 50th Anniversary of the Metropolis Algorithm at Los Alamos. To give proper credit to all the authors, it should be the Metropolis, Rosenbluth, Rosenbluth, Teller, and Teller (MRRTT) method, but most people just call it the "Metropolis method" (or algorithm).

I enjoyed meeting people at the conference, learning more about various techniques (look for future posts), and visiting New Mexico (my first time there).

One of the sessions was a history of the Metropolis method, and Marshall Rosenbluth gave a talk. It was good to hear from one of the original authors.

2003-05-12T00:02:00.000-05:00

Convergence Diagnostics

A while ago I expressed a desire for better automatic control of Monte Carlo simulations. The magic keyword seems to be "convergence diagnostics".

Statistics practitioners discovered that the Metropolis method (and other Markov Chain methods) were useful for integrals in Bayesian analysis. (handy overview)
I suppose being statistics researchers, they worried more about convergence more that physicists(?) Note there are two convergence worries - has the algorithm reached the long-time limit (stationary distribution), and has it converged to the desired distribution?

The most recent review of the various convergence diagnostics that I've found is by Brooks and Roberts (1999). Not all the diagnostics are applicable to Metropolis sampling (or Metropolis-Hastings, as generalized Metropolis is called).

The other interesting item to look at is "convergence acceleration" and see if any those techniques are useful for physics simulations.

2003-04-22T23:57:00.000-05:00

Randomness as a professional disease

In the May/June 2003 issue of Computing in Science and Engineering, Dietrich Stauffer has an article on Sociophyics Simulations (p. 71):

"Statistical physicists, of course, assume that these hierarchies arise from randomness (The illusion that everything is random is a professional disease - morbus Boltzmann - among such physicists, just as silicosis, or black lung, affects mine workers)"

2003-04-19T11:02:00.000-05:00

Generating trial moves cheaply

Today's article is "Monte Carlo simulations using sampling from an approximate potential" by Lev Gelb (JCP 118, 7747). He describes using an approxmiate potential to make a series of (cheap) moves and then accept or reject that series of moves based on the desired (correct, expensive) potential. This is similar to (but not quite the same as) Multilevel Monte Carlo in other work. The focus of multilevel sampling is early rejection of ultimately undesirable trial moves. The focus of Gelb's method seems to be creating uncorrelated trial moves to evaluate with the expensive potential.

I'm a little disappointed that the author didn't reference the similar multilevel Monte Carlo used in PIMC (RMP 67, 279) (see section V.F, p. 329) and VMC (JCP 113, 5123).

On page 7749, he mentions that an empirical potential could be used as the approximate potential and an ab initio potential for the correct potential. This has been done, sort of. (Also in Recent Advances in Quantum Monte Carlo Methods II (World Scientific 2002) or here .) I called the technique "pre-rejection". Consider using Gelb's method, but make only a single trial move using the approximate potential. If that trial move is rejected, then the energy change is clearly zero, and we know the move will be accepted without having to evaluate the expensive potential. This shortcut can be extended to more trial moves, but its value decreases dramatically as the probability that all the trial moves (with the approximate potential) will be rejected becomes very small.

Gelb applied his method to a Lennard-Jones fluid where the potential is cut and shifted. The approximate potential is the potential with a short cutoff and the correct potential has a larger cutoff.

Table I shows speedups and energies (and pressures and heat capacities) with error bars. There seems to a tradeoff - more speed up lowers the acceptance ratio and leads to larger error bars. This can be quantified by computing the efficiency 1/ (T * sigma^2) (or Speedup/sigma^2). Computing efficiencies for system 1 gives

System	Efficiency of Energy (x10⁴)
1(ref)	12
1(a)	14
1(b)	19
1(c)	14
1(d)	25

From the paper, it looks like a,b and c,d have the same parameters of relevance (rc and M) and that this method definitely does increase the efficiency.

2003-04-17T23:33:00.000-05:00

The Quest for Better Wavefunctions

One of the nice features of QMC is the freedom in the choice of wavefunction. The downside is it can be time consuming and unwieldy to optimize them in VMC.

In a recent article - Backflow Correlations for the Electron Gas and Metallic Hydrogen - Holzmann, et al. look at improving wavefunctions through adding analytic information.

The section with the Feynman-Kacs approach uses a cumulant expansion. In case you (read: me) need a reminder of the cumulant expansion, look at this online statistics text (chapter 8.4).

2003-04-09T23:42:00.000-05:00

Followup to "Testing Monte Carlo"

I wrote a program that uses Simpson's rule to compute the energy of two particles in a box interacting via a Lennard-Jones potential

The web page and program are here

It runs in a second or so on my 1 GHz P3. I expect 3 particles should be doable.

Edit 9/4/2004, fixed the link.

2003-04-07T20:28:00.000-05:00

Testing Monte Carlo

You've just written a brand new Monte Carlo code. It's producing vast amounts of data on your fast new computer. How do know if that random looking time series will average to a beautiful result and lead to deep insight, or is it simply random noise from the physics of some other universe?

There are several techniques we can use to verify that the codes we write are working correctly.

Unit testing
This is standard software engineering (or *should* be standard). This involves small drivers to test individual classes or routines. It's good for verifying small, easily understood pieces before integrating into a larger system and for regression testing (making sure the latest change didn't break something).
Solve simple systems with analytic or well-known answers.
The harmonic oscillator is popular. For QMC, H₂ is very well known, and very accurate answers are known from other techniques.
Compute internal quantities in multiple ways.
The best example in QMC is the local energy - compute the derivatives analytically and numerically. In PIMC there are multiple energy estimators.
Compute system values with another integration method.
It seems that computers are fast enough that we should be able to use grid-based methods to compute answers for small numbers (1-3) of particles.
For QMC, shutting off the electron-electron interaction simplifies the problem, and Mathematica can handle the resulting integrals (numerically, not symbolically)
Comparison with literature
Compare with experiment and answers obtained by other methods.

2003-03-28T23:44:00.000-06:00

Do you believe your error bars?

Computing uncertainty estimates from Monte Carlo data usually assumes the data is normal (or can be sufficiently reblocked to be normal, according to the central limit theorem). But what if we don't have quite that much data - what effect would the underlying distribution have on the uncertainty estimates (error bars)?

Mervlyn Moodley looks at this question in the paper The Lognormal Distribution and Quantum Monte Carlo Data.

The interesting figures seem to be figures 3 and 9. In Figure 3, we see the confidence intervals change as the data is reblocked (and hence comes closer to normal). Figure 9 is generated from a real simulation, and it seems there is very little difference between the error bars computed assuming normality and those taking the underlying lognormal distribution into account.

2003-03-17T22:14:00.000-06:00

Bilinear QMC

I'm reading a recent paper by Arias de Saavedra and Kalos titled "Bilinear diffusion quantum Monte Carlo methods" (PRE 67, 026708). They present an algorithm that uses a pair of walkers (the bilinear part) to sample the square of the ground state wavefunction (rather than first power of the wavefunction like DMC).

They claim it's useful for computing unbiased expectation values (see my post from Monday, Feb. 24 - although I'm not sure it would help with derivative operators) and energy differences.

The part I find most intriguing is the way they use importance sampling - they get reasonable results for hydrogen without it. And when they do use importance sampling, it's only to remove singularities due to cusps. I would guess that any successful scheme (ie, able to scale reasonably with system size) needs to use as an accurate a trial wavefunction as possible to guide the sampling.

The more I study the paper, the less I understand it. I guess the first step is to try reproduce their results for the 1-D harmonic oscillator.

2003-03-11T23:35:00.000-06:00

Steal this research idea

Automatic control of Monte Carlo simulations - I want to know how to do it, but I don't necessarily want to figure it out myself. :)

Determining equilibrium and phase transitions.
Typically the particles are started in a lattice or random configuration and let run a whlie to equilibrate. I want automatic detection of equilibrium, which is especially tricky near phase boundaries. Also near phase transitions, the system can suddenly change phases.
Adjustment of parameters for optimum efficiency
The key quantity affecting the efficiency is the correlation time, but it's very noisy. Is there some other quantity that we could use that would behave similarly, but with less noise? Since this will only affect the efficiency, approximations are quite okay. Like assuming some analytic form for the autocorrelation function with a small number of parameters and fitting to it?

Also note that SGA-like algorithms can be used for control as well as optimization. The decay parameter reaches a constant to follow the system, rather than continuing to decrease, as in optimization.
Better control of DMC population

Occasionally, I have trouble with the number of walkers increasing rapidly. In my codes, this means it exceeds some maximum value and stops the program. Sometimes the number of walkers jumps to a large but stable value, and the energy becomes far too low and clearly wrong. The hard part is that population control introduces a bias, and can't be too intrusive.

For these items, are there any concepts from control theory or signal processing that would be useful?

2003-03-04T23:02:00.000-06:00

Dynamic Monte Carlo

Checking the arxiv.org preprint server, there was an article titled An Introduction To Monte Carlo Simulations Of Surface Reactions by A.P.J. Jansen, which describes Dynamic Monte Carlo methods (Kinetic Monte Carlo is similar). One important input is transition probabilities. On page 21, "Estimates of the error made using DFT for such systems are at least about 10 kJ/mol."

Can QMC do better?

2003-02-28T00:33:00.000-06:00

VMC Optimization: SGA vs. fixed sample reweighting

In VMC, we parameterize the wave function and vary the parameters until the minimum energy (or variance) is reached. Fixed sample reweighting generates a number of configurations from the MC process, and then uses a least squares minimization algorithm and a reweighting estimate for the energy to minimize the variance of the energies.

The Stochastic Gradient Approximation (SGA) takes a noisy estimate of the gradient and moves in that direction ("downhill"). The amount it moves is scaled by a variable that decreases each iteration (usually 1/n). It is guaranteed to converge, but it may take a while.

I don't like the fixed sample reweighting (FSR) method on aesthetic grounds. SGA feels more in line with the "Monte Carlo spirit". FSR feels like it is taking very coarse steps: sample configurations for a while, stop and do a minimization, repeat these two steps for a few iterations. Each step is large, and results in significant improvement. The SGA, on the other hand takes a large number of small steps. Each individual step doesn't give much improvement (or any improvement at all), but the steps accumulate to give a good answer.

Aesthetic considerations aside, FSR requires no analytic gradients for efficient implementation. SGA needs analytic gradients for most of the parameters (you could do numerical gradients, but the time gets much worse)

FSR does variance minimization, and SGA does energy minimization - there may be some differences in the resulting quality of the wave functions (JCP 112, 4935)

2003-02-24T23:10:00.000-06:00

Derivative operators and forward walking in DMC

Using Diffusion Monte Carlo (DMC) to compute estimators other than the energy leads to a "mixed estimator"- an average that combines the trial wavefunction and the exact ground state (actually the fixed node ground state) - not what we really want. However, the number of children of each walker after a long time is proportional to the ratio of the fixed node ground state to the trial wavefunction. This leads to the forward walking procedure for obtaining correct averages over the fixed node ground state. See Casulleras and Boronat ( PRB 52, 3654) for examples, references, and implementation details.

Numerical values for derivative operators, such as the kinetic energy, are computed by applying the operator to the trial wavefunction. The forward walking procedure only corrects the distribution in the average, not the values of the estimator. (we really want the derivative applied to the exact wavefunction, not the trial wavefunction).

Computing the kinetic energy (and the pressure) by naively doing forward walking on the kinetic energy is incorrect. For the kinetic energy, we can use the total energy and the forward walked potential energy to obtain the kinetic energy.